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Abstract. The structure of poloidal and toroidal flows of trace impurities in the 
edge pedestal of tokamak plasmas is studied analytically and numerically. Parallel 
momentum balance is analysed upon retaining the following terms: poloidal and 
toroidal centrifugal forces (inertia), pressure force, electric force, and the friction force. 
It is shown that, when the poloidal flow is such to produce a properly defined Mach 
number of order unity somewhere on the flux surface, shock fronts can form. The 
shock fronts can modify the predicted asymmetry structures in both the flow and the 
density profile along the poloidal arc. Predictions of the theory are shown against 
experimental observations in the ASDEX Upgrade tokamak, showing good qualitative 
and quantitative agreement if the inertia term associated with the poloidal flow is 
retained. 
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1. Introduction 

Understanding the behaviour of impurities in the edge pedestal of tokamak H-modes 
is a fundamental requirement to be able to predict impurity penetration in the plasma 
core. The available theory, and it is referred here in particular to neoclassical theory, 
can very well predict flux-surface averaged transport fluxes and species poloidal flows 
and their relationship with the toroidal flow and the particle density, see for example 
PQ for a thorough review, and (2J |3] for more detailed calculations. 

It is well known that flows and density equilibration on a flux surface are regulated 
by collisional processes (mainly friction for trace impurities) and strong damping of 
acoustic modes perturbations. However, the exact shape of the density and flow profile 
on the flux surface is not provided by standard neoclassical theory. Previous works in 
this direction have been produced, assuming dominant impurity-on-bulk ions friction 
[I], or rigid toroidal rotation of the whole plasma [5j |6]- In Ref. [I] in particular 
the outcome of the predicted flow and density poloidal profiles on the radial transport 
is discussed. However, previous works have not considered the inertial term in the 
parallel momentum balance. As it turns out and it will be shown later on, that term is 
necessary when dealing with flows that can become super-sonic along the poloidal arc. 
In fact, standard neoclassical ordering breaks down for sonic flows, as in that case the 
perturbation to a 'global' Maxwellian is of the same order of the equilibrium function. 

Previous works on the impact of sonic poloidal flows on the bulk species behavior 
and on the magnetic equilibrium have been published [8], however here the focus 
is on the impurities, while bulk species are assumed sub-sonic in flow speed. In this 
work it is shown, through both analytical and numerical calculations, how the poloidal 
structure of the flow and density of a trace impurity emerges by the interplay of inertia, 
friction, and the pressure and electric drive. Since the magnetic field acts effectively as 
a nozzle, it can be easily seen that static shock fronts can form as the poloidal Mach 
number, properly defined, goes through one. Peculiar of the toroidal plasma is that the 
nozzle is periodic and radially coupled, meaning that shock fronts at different poloidal 
positions appear at the different radial locations, implying that radial coupling and a 
three-dimensional flow structure can form around the shock fronts [9]. 

The present work has been inspired following charge exchange measurements 
[T0| ITT] and observations of poloidal asymmetries of impurity flows in the edge pedestal 
[T2l [TBI ITT] in the ASDEX Upgrade (AUG) tokamak. Observations on other machines 
have been reported as well, see e.g. [T5] . 

The manuscript is organized as follows: in section 2 the main assumptions and 
equations are introduced. In section 3 the main results are derived analytically. In 
section 4 numerical examples are computed and qualitatively compared to experimental 
observations. Section 5 draws the conclusions. 
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2. Main assumptions and equations 

2.1. Some definitions and assumptions 

The background magnetic field is defined by its toroidal component = FR (F being 
the diamagnetic flux function and R the local major radius) and its poloidal component 
Bg = | Vi/j\/R, where is the poloidal magnetic flux and R is the local major radius. An 
appropriate coordinate system is defined as the triplet (ip , 9 , <p) , with 9,(f> respectively 
the poloidal and toroidal angle (with 9 increasing upwards from the outboard side and (f> 
increasing clockwise when the torus is viewed from above). The actual signs of the fields 
are also introduced: a^, <jq respectively for the toroidal and the poloidal field. As such 
in the following it assumed strictly that B^^B e > 0, while the signs are explicited in 
the equations through as. A reference location on the poloidal angle is chosen as 9 = 0, 
the outboard mid-plane position, and will be identified in the following with the pedix 
'L', while the inboard mid-plane position on the high field side will be identified with 
the pedix 'H' . The usual flux-surface-average (FSA) is defined as (/) = d/dV J fdV, 
where V is the plasma volume enclosed by the local flux surface. It is reminded a basic 
property of the FSA: (B • V/) = 0. 

The main ion species Mach number is defined as M a>i = |Uj|/^ h , with U; the 
species flow and v [ h = y/Ti/Mi the species thermal speed. The impurity Mach number 
is defined as M^ z = |Uz|/f t z h , with XJ Z the impurity velocity and t> t z h = ^Ti/M z the 
impurity thermal speed. In the following it is assumed that bulk ions are strongly 
subsonic, M &) j <C 1, while there is no limit to the speed at which the trace impurities 
can rotate. 

It is now clarified what is the definition of 'trace' as used in the following for the 
impurities: it is assumed that Zn z <C Ui, where n z is the impurity density and n; the 
main ions density. This means that the background geometry and electrostatic potential 
are determined by the bulk species only. On the other hand, to be able to treat the 
impurities in a fluid context, it is assumed that the parallel Knudsen number, defined 
as K n ^ = vf h /(Rqu zz ) < 1, with v zz the impurity self-collisional frequency and Rq 
the parallel connection length. Also, the Knudsen number associated to the impurity- 
ion collisions is given by if n ,||,z-i — v^ h / \Rqv Z i) < 1. Both Knudsen numbers being 
below unity, means that the impurity distribution function can be represented as a local 
Maxwellian, which describes an isotropic pressure tensor, and for which the parallel fluid 
flow will be influenced mostly by friction with the bulk ions. 

Finally, it is assumed that the ions species temperature is a flux function, T; L = 
Ti(ifj), while the impurity density can acquire poloidal dependence. The main ion density 
is assumed to be a flux function as well, since it is considered here a situation in which 
M M < 1. 

In the remainder, the pedix Z characterizing the impurity species is implicitly 
assumed for density, temperature, and velocities. It will be explicited when needed (e.g. 
for main ions or electrons). 
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2.2. Particle continuity equation for the impurity 

The particle continuity equation is written (for impurity species) as: 

V • (nU) = (1) 

where the d/dt term is neglected (it is of the order of the cross-field flow). The flow is 
predominantly on a flux surface, i.e. |U ■ e^| <C |U|. In this case equation (JTJ) reduces 
to the well known poloidal dependence of the poloidal flow component U$: 

where K is a yet to be determined flux function. Working on the general velocity 
expression U = U\\h + XJ± = U^e^ + Uee p , one finds the toroidal and parallel velocity 
components expressed as: 

= a^agK— + a e RVt 

n 

B Ftt 

U\\ = a e K h oea^— - (3) 

n B 

with Q = —d&o/dijj the rigid rotator speed. Note that it is assumed here an 'ideal-like' 
Ohm's law E + U x B = 0, with E = — V$ - The diamagnetic term oc VP/{Zn) 
is neglected at this stage. As such the electric field is computed from the known flow 
velocities at one specific poloidal position. The role of the diamagnetic flow will not be 
investigated further and is left for future work. The zero-order electrostatic potential 
3>o(V0 is a A ux function. As such also Q is a flux function. A relation between K and 
the FSA parallel velocity can also be found: 

(BU\\) Ftl 

K = ae WM " a "WM (4) 

Note that a pure toroidal flow can be sustained by choosing (BU\\) = aea^Ffl, however a 
poloidal flow is always accompanied by a toroidal flow with in-out asymmetries (Pfirsch- 
Schluter flow). 



2.3. Parallel force balance for impurities 
The parallel force balance equation is simply: 

MnB ■ [(U ■ V) U] = -TBVfi - ZnBV\\® + BT\\ (5) 

where the self-viscous term and the inductive electric force are neglected, the first for 
the trace assumption, the second since most H-modes work at low edge loop voltage. J-u 
is the parallel component of the friction force, which is considered between the impurity 
and the main ions only: 

JT,, = -Mnv Zi {U\\ - U u ) (6) 

where vz\ is the impurity-to-main ions collisional frequency. The poloidally varying 
potential $ is assumed as a small perturbation with respect to the dominant one $o- 
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The inertial term, i.e. the first term on the left-hand-side of equation (jSJ), can be 

fU 2 \ 

manipulated as (U ■ V) U = V I — -Ux(VxU). Upon substituting formulas pj), 



x 2 r 

after some algebra, the following equation is obtained: 

MB ' V (l^ ~ = ~ TBV W hgn ~ ZB ^\\® - Mv n B{U\\ - [7i.ii) (7) 

As a remark, notice that the collisional frequency, being z/zi oc n; is a flux function. 

This fact allows to compute the averaged parallel velocity from the FSA of equation (0) 

simply as: (BU\\) = ay = (BUi\\), leading to K = 6 V °\ . The main ions 

\B l jn) 

averaged parallel flow term ay is computed using information from the flow profiles 
obtained at one specific poloidal location of reference. 

Equation (J7J) is finally recast in Bernoulli-like form (v t ^ = y^T/M): 



B V 



K 2 B 2 „ 2 2 
R 2 Vl 2 + 2v 2 



7? 2 




-2^5(^,1-^,11) (8) 



2.4- Pfirsch-Schlueter $ps and centrifugal potential $ w 

To evaluate the source $, to enter in equation (JE1), it is first followed Ref. [IE], which 
deals with the Pfirsch-Schliiter potential, leading to: 

$ps = a^FP' f^(l- (9) 



L 



Be \ g 3 R 2 , 

Since P' is negative, and at the LFS both Bg is smaller and (#3 — 1/R 2 ) is positive, it 
turns out that $ is negative at the HFS equatorial mid-plane. 

In addition, one can also include the centrifugal potential |6j, which is given by 

\ uj 2 R\ (R 2 

where u is the toroidal angular rotation frequency of the main ions. The resulting 
potential perturbation is thus $ = $ps + "3V 



= 155- - 1 ( 10 ) 



2.5. Phenomenological diffusion 

It is shown in the next section that equation (jSJ) can develop shock fronts [S]. Since 
along the radial direction, these shocks appear at different poloidal locations (basically 
the shock front poloidal location is determined by the density and the collisionality) , 
it is expected that a smooth structure finally is sustained [9]. Note that gyrokinetic 
turbulence is not a candidate here since turbulence structures are of the order of the 
Larmor radius. On the other hand, parallel transport due to free streaming can also 
smooth the shock front. 



Impurity flow structure in the tokamak edge pedestal 



6 



2.6. Consistency with the gyrokinetic equilibrium 

It is now shown that equation (jHJ) can be derived from the equilibrium gyrokinetic 
equation, demonstrating its consistency. The gyrokinetic equation for the impurity 
equilibrium distribution function F is: 

V- [(v\\b + v E )BF] (y\\BF) =C (11) 

where ve = Ex B/B 2 , E = — V$o, and C is the collision operator. All terms inversely 

proportional to the charge Z have been neglected (i.e. the curvature and centrifugal 

ve • V.B 

drifts). The parallel acceleration is given by: vu = — /ib- V-B— Ze/Mh- V$+W|| — h 

B 

ve-[(ve • V) b], where it is important to retain both terms involving ve- The distribution 
function F is assumed as an isotropic local Maxwellian F oc ne -0 ' 5 ^ 11- ^ 11 ^ and for 
reference it is set (mass, temperature) M = T = 1. 

First integration of the equation in velocity space J -dfxdvu leads to the continuity 
equation V ■ (nUnb + nv E ) = 0. Then, integrating with J -vnd/idvn leads to the parallel 
force balance equation: 

V • (nU*b + nb + nU\\V E ) + 

^b • VB + Znb ■ V<f> - nU\\ WE - nv E • [(v E • V) b] = Ru (12) 
B B 

where R\\ is the parallel friction term. After some algebra it is finally found that 

equation (fl2i) coincides with equation (jHJ), where the relevant identity ve ■ [(ve ■ V) b] = 
~R 2 n 2 Bl 



b- V 



has been derived. 



2 B 2 

3. Poloidal flow structure and density distribution of the trace impurity 

Equation (jSJ) is now studied in detail by first putting it into dimensionless form as: 

A [ a f - r 2 + (3 ( 7 - log/)] = v x + u 2 f + (13) 

where the unknown / = (B 2 /^^ 1 B 2 /n, and the other parameters are defined as: 
a = K(B 2 /n)/(2B 2 v th ), r = R 2 Q 2 /(2K(B 2 /n)v th ), (3 = v th /(K(B 2 /n)), 7 = logS 2 + 
Z$/r, v x = 5(a e auB 2 - a^FtlB 2 ) /(B 2 K(B 2 /n}), v 2 = -5{a e au - a^FQ) / (K (B 2 /n)), 
5 = RlVzi/ (BoVth). The choice of the independent variable / is such that (/) = 1 is 
an automatic constraint of the problem, additionally to / being periodic in x after one 
poloidal turn. The coordinate x is defined such that dx = dd j '(i?£,|V0|). The physical 
meaning of the different terms is respectively: (left-hand-side) poloidal centrifugal 
energy, toroidal centrifugal energy, nozzle striction effect from magnetic field and electric 
potential energy, pressure drive, (right-hand-side) friction breaking from toroidal and 
poloidal rotation of the main ions, and finally the phenomenological diffusion term. 

The mathematical structure of equation ( {TBI is that of a viscous Burger equation 
with external forcing, while the case /i = is equivalent to the inviscid Burger equation 
with additional linear contributions. The latter, inviscid limit, can be rewritten very 
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compactly as (2a f 2 — /3) — — V af 3 + bf 2 + cf = d. The latter form shows clearly that 

a singularity in the equation appears when / = y//3/(2a), except if, at the same time, 
for said value of /, af 3 + bf 2 + cf — d = 0. In principle such a location, where |^ is 
finite despite its coefficient vanishing, exists along the poloidal arc. On the other hand, 
because the system is periodic, there has to be another location at which both conditions 
will not be satisfied. This results in a shock front appearing, producing a finite jump in 
/, i.e. in density, poloidal and toroidal flow. This phenomenon will always appear as 
soon as a location where / = a//?/ (2a) is present. Note that said condition is physically 
the fact that the Mach number associated with the quantity UqB/Bq is unity. In usual 
H-mode edge plasmas such Mach numbers can be easily reached for the impurities, 
which are not subject to poloidal flow damping when being in trace concentration. 

As example of analytical limits of equation ( fl3|) . consider the cases > a and 
low collisionality, then one obtains / ~ B 2 , which results in the density n being a flux 
function. Another case is when the toroidal centrifugal term r is dominant, leading to 
the well known outboard impurity accumulation [6]. The limit of infinite collisionality 
as studied in [3], leads to / = 1, i.e. n ~ B 2 , which provides inboard impurity 
accumulation. The limit found when, neglecting collisions, a>)5, leads again to / ~ 1, 
i.e. n pd B 2 . In this respect a super-sonic collisionless flow and a strongly collisional 
flow are indistinguishable from the point of view of the density poloidal distribution. On 
the other hand, when a ~ ft, with or without collisions, a peculiar structure appears, 
i.e. the shock front, which has observable characteristics. Realistically most edge El- 
mo de plasmas are in this regime of mixed transonic-collisional flow where all terms in 
equation f Tl3|) are important. 

4. Numerical applications and comparison with experimental observations 

Equation (TT3]) is solved numerically by employing a numerical scheme based on 
'convergent time evolution', i.e. adding a time derivative term on the left-hand-side 
and make it relax through numerical diffusion, and spatially by using upwind Euler 
scheme which is automatically stable for large flow velocities, however it has intrinsic 
numerical diffusion. The upwind Euler scheme nicely resolves the static shock front and 
automatically stabilizes the artificial 'sound' waves generated by the convective non- 
linearity. Iterations are also employed to solve for the unknown term (B 2 /n). The time 
step is chosen to reach steady-state and the spatial grid is sufficiently large to resolve the 
poloidal structures and the diffusion layer (of the order of yfii). The numerical scheme 
has been validated by benchmarking against analytical expressions of the equation in 
different limits and by checking conservation of (/) = 1. 

The input profiles are taken from charge exchange measurements of both LFS and 
HFS poloidal and toroidal rotation in ASDEX Upgrade [11] , for boron B 5+ . Also the 
LFS impurity density and temperature are measured and used as input. The values 
of au and Q are computed from the measured flow velocities, while the deuterium 
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density, relevant for the impurity-ion collisionality, is taken as having the same impurity 
density profile, but rescaled such that the pedestal top value matches the one from 
the experimentally measured electron density. The magnetic fields are taken from the 
equilibrium reconstruction code CLISTE [17]. Note that in the standard magnetic 
configuration of AUG — +1, a$ — — 1. Profiles of the input quantities can be found in 
figure [Ua,b,c). In figure [Qa), a modified poloidal velocity is also used (dot-dashed line), 
which has constant velocity from the position of the maximum (p^ « 0.99) up to the 
separatrix. The reason for using this modified velocity profile is to somewhat consider 
eventual uncertainties in separatrix position with respect to the rotation measurements 
and in experimental uncertainties. 

In figures 2, 3, and 4 the results of the calculations are shown, where in subplots (a) 
the inertia term proportional to a and r, as well as diffusion, are set to zero, wheareas 
in subplots (b) the equation is solved in its entirety. Here the viscosity parameter is 
chosen of the order of /i ~ 0.01, which does not disturb the resulting profiles but helps in 
smoothing the shock front. At present the exact value of fi can not be determined from 
principles, thus requiring in future to perform global simulations as for example proposed 
in Ref. [20J. Notice that in figure EJ^c), the modified poloidal velocity of figure (U^a) has 
been employed. In figure [3](c) the poloidal distribution of the resulting toroidal flow is 
also shown, for the case with all terms retained in the numerical calculation. 

From the results displayed in the figures, several points can be clarified: 

(1) while the friction force alone can give a reasonable agreement in the poloidal flow 
asymmetry (not shown), it fails in reproducing the toroidal flow asymmetry, as can be 
seen in figure Et a); 

(2) the "poloidal" centrifugal force, at Mach of order unity, provides that additional 
poloidal structure that leads to matching both the poloidal and toroidal flow radial 
variations on LFS and HFS, as can be seen from figure Etb,c). Notice also that the LFS 
value of the poloidal flow regulates the behavior of the toroidal flow asymmetry, as when 
the former quantity drops, also the latter behave differently, as comparing figure EJ^b) 
and figure Etc); 

(3) there is a strong indication that at the top of the plasma the toroidal flow has 
reversed direction with respect to both LFS and HFS, as can be seen in figure EJ^c). 
This is a testable prediction of the calculations presented here; 

(4) smoothing of the shock structure is essential, as seen in comparing figure @|a) and 
figure |4jb). 

5. Conclusions 

The poloidal structure of flows and density of trace impurities is studied by employing 
parallel momentum balance and retaining the non-linear convective (or inertial) terms 
that produce additional effects with respect to previous results [HE]. In particular 
the non-linear term arising from the poloidal flow is necessary to properly introduce 
shock fronts and poloidal shock structures that drastically change the resulting impurity 
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density profiles when the associated flow goes from sub to super-sonic along the poloidal 
arc. 

The findings reported here show that it should be rather generic, in the H-mode 
edge, for these structures to exist since trace impurities can be strongly accelerated in 
the poloidal direction by the background electric field even in the absence of average 
parallel rotation of the main ions. A remark has to be made with respect to the type 
of impurity, since the effects presented here should be valid for low and moderate- Z 
impurities, while for high-Z impurities the role of the poloidally asymmetric potential 
and of the toroidal centrifugal force could become important and modify the outcome 

[MEE]. 

Comparison of the discussed theory against experimental observations in ASDEX 
Upgrade is showing that all terms are important to reproduce the observed flow 
structures, at least on the LFS and HFS. Measurements on the top of the plasma 
should provide a key answer since the shock front will deeply influence the LFS-top- 
HFS profile. 

Additional effects that can play a role in the pedestal and are not considered here 
are: diamagnetic flow, local sources, orbit loss, global effects, application of localized 
radio-frequency heating. It is expected that the orbit loss mechanism in particular 
would be important since the pedestal region has a width which is usually of the order 
of the impurity poloidal Larmor radius, as well as global effects [20J. 
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Figure 1. Radial profiles with error bars, on the coordinate, of 
experimentally measured (a) LFS poloidal rotation Uq (solid), HFS poloidal 
rotation (dashed), and (b) LFS (dashed) and HFS (dashed) toroidal rotation 
U^; (c) LFS impurity density. In subplot (a), the dot-dashed line is the 
poloidal velocity, modified such that the value stays constant from the position 
of maximum value up to the separatrix. This example is taken from ELMy 
H-mode AUG discharge #28093. 
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Figure 2. (a) Comparison of radial profiles of LFS and HFS toroidal rotation for 
the case with a = fi = 0; (b) same comparison but including all terms; (c) same 
comparison but here the input poloidal rotation profile has a constant value 
from the position of the maximum velocity (p^ ~ 0.99) up to the separatrix. 
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Figure 3. (a) Poloidal variation of impurity density for different radial locations 
(in legend) for the case with a = [/, = 0; (b) same plot but with all terms 
included; (c) Poloidal distribution of the toroidal rotation profile for different 
radial locations (in legend), case with all terms included. 




Figure 4. (a) Poloidal variation of Mach number when both friction and 
viscosity are neglected for three radial locations (in legend); (b) same plot 
but with all terms included. 



